System, method and computer program product for modeling at least one section of a curve

ABSTRACT

A system, method and computer program product are provided for modeling at least one section of a curve. Each section can be modeled by initially providing a pair of positions (f i , f i+1 ) of the section of the curve including associated directions (d i ,d i+1 ) and curvatures (κ i ,κ i+1 ). Then, points b 0 ,b 1 ,b 2 ,b 3  and b 4  can be identified based upon the pair of positions (f i ,f i+1 ) and associated directions (d i ,d i+1 ) and curvatures (κ i ,κ i+1 ). Thereafter, a quartic interpolant p(t) can be determined over an interval (i≦t≦i+ 1 ) based upon points b 0 ,b 1 ,b 2 ,b 3  and b 4  to thereby model the section of the curve. The quartic interpolant can be determined such that the interpolant p(t) has a position, direction and curvature equal to f i , d i  and κ i , respectively, at t=i, and the interpolant p(t) has a position, direction and curvature equal to f i+1 , d i+1 , and κ i+1 , respectively, at t=i+1.

FIELD OF THE INVENTION

[0001] The present invention relates generally to systems, methods and computer program products for modeling curves and, more particularly, relates to systems, methods and computer program products for modeling at least one section of a curve according to a quartic Hermite interpolation technique.

BACKGROUND OF THE INVENTION

[0002] In a number of different disciplines, such as computer-aided geometric design (CAGD), it is necessary to model, or approximate, parametric curves from sparse data. Typically, such modeling techniques attempt to model planar curves with high accuracy, while maintaining curvature continuity. In this regard, many conventional techniques model curves based upon interpolation algorithms that are based upon the notion of the curve being twice continuously differentiable with respect to arclength and capable of being constructed locally.

[0003] A planar curve f is said to be twice continuously differentiable if the unit tangent vector f*:=f′/|f′| and the signed curvature f**:=f′×f″/|f′|³ are continuous. Let f_(i), d_(i) and κ_(i) be prescribed positions, directions and curvature values, respectively. For example, these data might be thought to have come from a twice continuously differentiable planar curve as:

f _(i) :=f(t _(i)), d _(i) :=f*(t _(i)) and κ_(i) :=f**(t ₁)  (1)

[0004] From such data, then, one such interpolation technique, known as Hermite interpolation, can determine a geometric Hermite interpolant p that satisfies:

p(i)=f_(i) , p*(i)=d _(i) and p**(i)=κ_(i)  (2)

[0005] where the components of p are polynomials on each parameter interval [i, i+1]. Also, to preserve the convexity of the data, that is, to satisfy the following assumption:

sign(κ_(i))=sign(d _(i)×(f _(i+1) −f ₁))=sign(f _(i) −f _(i−1))×d _(i)), for all i  (3)

[0006] the curvature of the interpolant typically has the same sign throughout.

[0007] As will be appreciated by those skilled in the art, many modeling, or interpolation techniques, operate to model a curve by separately interpolating different sections of the curve, where each section is defined by a set of end points. In this regard, because p can be constructed interval-by-interval, any such piecewise constructed interpolant will necessarily be continuously differentiable. Therefore, the interpolant p can be described as being constructed over the interval [0,1] without any loss of generality. According to one Hermite interpolation technique, developed by de Boor, Hollig and Sabin (referred to herein as the BHS technique), a cubic polynomial is used on each interval. Turning attention to the interval [0,1], then, the interpolant p can be written as: $\begin{matrix} {{{p(t)} = {\sum\limits_{j = 0}^{3}{b_{j}{B_{j}(t)}}}},} & \left( {0 \leq t \leq 1} \right) \end{matrix}$ ${{with}\quad {B_{j}(t)}}:={\begin{pmatrix} 3 \\ j \end{pmatrix}{{t^{j}\left( {1 - t} \right)}^{({3 - j})}.}}$

[0008] Interpolation of position at the given knots implies that

[0009] b₀=f₀ and b₃=f₁.

[0010] The remaining interpolation conditions can be simplified by introducing $\overset{\sim}{b}:=\frac{{\left( {f_{1} \times d_{1}} \right)d_{0}} + {\left( {d_{0} \times f_{0}} \right)d_{1}}}{d_{0} \times d_{1}}$

[0011] Geometrically, {tilde over (b)} is the intersection of the line through f₀ parallel to d₀, and the line through f₁ parallel to d₁, as shown in FIG. 1. The tangent interpolation condition is then defined as:

b ₁=(1−ρ₀)b ₀+ρ₀ {tilde over (b)} and b ₂=(1−ρ₁)b ₃+ρ₁ {tilde over (b)}

[0012] for some 0<ρ₀, ρ₁≦1, where ρ₀ and ρ₁ represent weights in writing b₁ and b₂ as convex combinations of b₀, b₃ and {tilde over (b)}. Moreover, the curvature interpolation condition is given by: $\begin{matrix} {\kappa_{0} = \frac{2d_{0} \times \left( {b_{2} - b_{1}} \right)}{3{{b_{1} - b_{0}}}^{2}}} & {and} & {\kappa_{1} = {\frac{2\left( {b_{2} - b_{1}} \right) \times d_{1}}{3{{b_{3} - b_{2}}}^{2}}.}} \end{matrix}$

[0013] Substituting b₂−b₁=b₂−{tilde over (b)}+{tilde over (b)}−b₁=(1−ρ₁)(b₃−{tilde over (b)})+(1−ρ₀)({tilde over (b)}−b₀) into the last equations for the curvature interpolation condition, κ₀ and κ₁ can be determined as follows: $\begin{matrix} \begin{matrix} {\kappa_{0} = \frac{2\left( {1 - \rho_{1}} \right)d_{0} \times \left( {b_{3} - \overset{\sim}{b}} \right)}{3\quad \rho_{0}^{2}{{\overset{\sim}{b} - b_{0}}}^{2}}} & {and} & {\kappa_{1} = \frac{2\left( {1 - \rho_{0}} \right)\left( {\overset{\sim}{b} - b_{0}} \right) \times d_{1}}{3\quad \rho_{1}^{2}{{b_{3} - \overset{\sim}{b}}}^{2}}} \end{matrix} & (4) \end{matrix}$

[0014] Defining an R₀ and R₁ as follows: $\begin{matrix} {R_{0}:=\frac{3\quad \kappa_{0}{{\overset{\sim}{b} - b_{0}}}^{2}}{2d_{0} \times \left( {b_{3} - \overset{\sim}{b}} \right)}} & {and} & {R_{1}:=\frac{3\quad \kappa_{1}{{b_{3} - \overset{\sim}{b}}}^{2}}{2\left( {\overset{\sim}{b} - b_{0}} \right) \times d_{1}}} \end{matrix}$

[0015] and solving equation (4) for ρ₀ and ρ₁ gives:

ρ₀=1−R ₁ρ₁ ² and ρ₁=1−R ₀ρ₀ ².

[0016] As will be appreciated by those skilled in the art, the BHS technique may have 0, 1, 2 or 3 solutions with 0<ρ₀, ρ₁≦1 and a sufficient condition for solvability is that (1−R₀)(1−R₁)≧0. With respect to techniques such as the BHS technique, it has been suggested that f_(i) and d_(i) be prescribed first, with appropriate values for κ_(i) selected to satisfy the foregoing condition. However, such a technique of determining f_(i), d_(i) and κ_(i) is not always feasible since, for example, there may be predetermined values that κ_(i) must have. Also, such a condition, which in terms of the data can be represented as:

[0017] (2(d₀×(f₁−f₀))(d₀×d₁)²−3κ₀((f₁−f₀)×d₁)²).

(2((f₁−f₀)×d₁)(d₀×d₁)²−3κ₁(d₀×(f₁−f₀))²)≧0,

[0018] puts an abstruse restriction on the data. Further, ρ₀ and ρ₁ cannot be chosen in any manner that will make them continuous functions of the data for all admissible data sets. For example, as shown in FIG. 2, when the data changes in such a way that R₁ remains constant at 0.875 while R₀ increases to 0.875 (say from 0), then ρ₀ =0.1674 and ρ₁=0.9755 when R₀=R₁=0.875. However, if the same data configuration is arrived at in such a way that R₀ remains constant at 0.875 while R₁ increases to 0.875 (say from 0), then ρ₀=0.9755 and ρ₁=0.1674 when R₀=R₁=0.875.

SUMMARY OF THE INVENTION

[0019] In light of the foregoing background, the present invention provides an improved system, method and computer program product for modeling at least one section of a curve. The system, method and computer program product of embodiments of the present invention are capable of modeling each section of a curve f based upon data including a pair of positions (f_(i), f_(i+1)) of the section of the curve, wherein the pair of positions includes associated directions (d_(i), d_(i+1)) and associated curvatures (κ_(i), κ_(i+1)). Advantageously, embodiments of the present invention are capable of modeling each section to exactly match the positions, directions and curvatures of data provided for the respective section. The system, method and computer program product of embodiments of the present invention are also capable of modeling each section of the curve to thereby preserve the shape properties of the curve f, such as by maintaining the convexity in the modeled section.

[0020] As indicated above, the system, method and computer program product of embodiments of the present invention are capable of modeling each section of a curve. Advantageously, then, changes in the data only affect the curve locally. That is to say, altering one data position (e.g., f₁) only alters the modeled curve up to adjacent data point(s) (e.g., f_(i+1)). As such, a point in one section can be altered without affecting the modeled curve in subsequent sections. In addition, the system, method and computer program product of embodiments of the present invention are capable of modeling the curve such that a change in the model is directly proportional to a change in the data. In other terms, the model depends continuously on the data such that small changes in the data (position, direction and/or curvature) will produce only small changes in the resulting model.

[0021] The system, method and computer program product of embodiments of the present invention are further capable of modeling each section of the curve with sixth-order accuracy. As will be appreciated, due to the finite number of data points that may be provided, no modeled section of the curve can exactly reproduce the respective section of the curve f. In this regard, the accuracy of a model can be defined based upon how much closer the modeled section comes to the respective section of the curve as the number of provided positions from the curve f increase. As such, by modeling each section of the curve with sixth-order accuracy, as the number of provided positions from a section of the curve increase by a factor of g, the modeled section of the curve will move g⁶ closer to the section of the curve.

[0022] According to one aspect of the present invention, a method is provided for modeling at least one section of a curve. According to one embodiment, a section is modeled by initially providing a pair of positions (f_(i), f_(i+1)) of the section of the curve, where the pair of positions includes associated directions (d_(i), d_(i+1)) and associated curvatures (κ_(i), κ_(i+1)). Then, points b₀, b₁, b₂, b₃ and b₄ are identified based upon the pair of positions (f_(i), f_(i+1)) and associated directions (d_(i), d_(i+1)) and curvatures (κ_(i), κ_(i+1)). More particularly, points b₀, b₁, b₂, b₃ and b₄ can be identified by initially defining a control point {tilde over (b)} based upon the pair of positions (f_(i), f_(i+1)) and the associated directions (d_(i), d_(i+1)). The control point can be defined as a point at the intersection of a line through point f₀ parallel to direction d₀, and a line through point f₁ and parallel to direction d₁. In other terms, the control point can be defined according to the following: $\overset{\sim}{b}:={\frac{{\left( {f_{1} \times d_{1}} \right)d_{0}} + {\left( {d_{0} \times f_{0}} \right)d_{1}}}{d_{0} \times d_{1}}.}$

[0023] After defining the control point {tilde over (b)}, points b₁ and b₃, can be identified. In this regard, point b₁ can be identified on a segment {overscore (b₀{tilde over (b)})} and point b₃ can be identified on a segment {overscore (b₄{tilde over (b)})}, where b₀ equals f_(i) and b₄ equals f_(i+1). For example, points b₁ and b₃ can be identified by defining weights (ρ₀, ρ₁) such that 0<ρ₀, ρ₁≦1. Thereafter, points b₁ and b₃ can be identified according to the following:

b ₁=(1−ρ₀)b ₀+ρ₀ {tilde over (b)} and b ₃=(1−ρ₁)b ₄+ρ₁ {tilde over (b)}.

[0024] Also after defining the control point {tilde over (b)}, point b₂ can be identified such that b₂ is disposed within a region bounded by a triangular shape defined by points {tilde over (b)}, b₁ and b₂, where point b₂ is identified based upon the curvatures (κ_(i), κ_(i+1)). Advantageously, by so identifying point b₂, the convexity of the modeled section can be preserved. For example, point b₂ can be identified by defining weights (α₀, α₁) based upon the curvatures (κ_(i), κ_(i+1)), where the weights (α₀, α₁) are defined such that α₀, α₁≧0 and α₀+α₁≦1. Thereafter, point b₂ can be identified according to the following:

b ₂=α₀ b ₁+α₁ b ₃+(1−α₀−α₁){tilde over (b)}.

[0025] In one embodiment, weights (α₀, α₁) can be defined by defining weights (ρ₀, ρ₁) and thereafter defining weights (α₀, α₁) according to the following: $\begin{matrix} {\alpha_{0} = \frac{R_{1}\rho_{1}^{2}}{1 - \rho_{0}}} & {and} & {{\alpha_{1} = \frac{R_{0}\rho_{0}^{2}}{1 - \rho_{1}}},} \end{matrix}$

[0026] where $\begin{matrix} {R_{0}:=\frac{4\quad \kappa_{i}{{\overset{\sim}{b} - b_{0}}}^{2}}{3d_{i} \times \left( {b_{4} - \overset{\sim}{b}} \right)}} & {and} & {R_{1}:=\frac{4\quad \kappa_{i + 1}{{b_{4} - \overset{\sim}{b}}}^{2}}{3\left( {\overset{\sim}{b} - b_{0}} \right) \times d_{i + 1}}} \end{matrix}.$

[0027] After identifying points b₀, b₁, b₂, b₃ and b₄, a quartic interpolant p(t) is determined over an interval (i≦t≦i+1) based upon points b₀, b₁, b₂, b₃ and b₄ to thereby model the section of the curve. More particularly, the interpolant p(t) can be determined according to the following: ${{p(t)} = {\sum\limits_{j = 0}^{4}{b_{j}{B_{j}(t)}}}},$

[0028] where ${B_{j}(t)}:={\begin{pmatrix} 4 \\ j \end{pmatrix}{{t^{j}\left( {1 - t} \right)}^{({4 - j})}.}}$

[0029] Advantageously, and illustrating that embodiments of the present invention are capable of modeling each section to exactly match the positions, directions and curvatures for the respective section, the interpolant p(t) has a position, direction and curvature equal to f_(i), d_(i) and κ_(i), respectively, at t=i, and the interpolant p(t) has a position, direction and curvature equal to f_(i+1), d_(i+1), and κ_(i+1), respectively, at t=i+1.

[0030] According to other aspects of the present invention, a system and computer program product are provided for modeling at least one section of a curve f.

BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Having thus described the invention in general terms, reference will now be made to the accompanying drawings, which are not necessarily drawn to scale, and wherein:

[0032]FIG. 1 illustrates a Bezier polygon with a corresponding cubic curve segment;

[0033]FIG. 2 illustrates a situation in which, ρ₀ and ρ₁ cannot be chosen to be continuous functions of data of planar curve f for all admissible data sets, according to one conventional interpolation technique;

[0034]FIG. 3 is a flow chart illustrating various steps in a method of modeling each of at least one section of a curve according to one embodiment of the present invention;

[0035]FIG. 4 illustrates a Bezier polygon with a corresponding quartic curve segment according to one embodiment of the present invention;

[0036]FIG. 5 illustrates the geometry of data from a circular arc according to one embodiment of the present invention;

[0037]FIG. 6 illustrates a plot of ρ_(opt)(θ)/ρ_(cnt)(θ) for −π<θ<π;

[0038]FIG. 7A illustrates a set of positions including associated directions and curvatures according to one embodiment of the present invention in the context of designing a dimensionless airfoil;

[0039]FIG. 7B illustrates a curve modeled based upon the set of positions from FIG. 7A according to one embodiment of the present invention; and

[0040]FIG. 8 is a schematic block diagram of the system of one embodiment of the present invention embodied by a computer.

DETAILED DESCRIPTION OF THE INVENTION

[0041] The present invention now will be described more fully hereinafter with reference to the accompanying drawings, in which preferred embodiments of the invention are shown. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein; rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art. Like numbers refer to like elements throughout.

[0042] As indicated above in the background section, a planar curve f is said to be twice continuously differentiable if the unit tangent vector f*:=f′/|f′| and the signed curvature f**:=f′×f″/|f′|³ are continuous. Again, let f_(i), d_(i) and κ_(i) be prescribed positions, directions and curvature values, respectively. For example, these data might be thought to have come from a twice continuously differentiable planar curve as:

f ₁ :=f(t _(i)), d _(i) :=f*(t ₁) and κ_(i) :=f**(t ₁)  (1)

[0043] According to one embodiment of the present invention, then, a system, method and computer program product are provided for modeling at least one section of a curve, where each section is modeled by determining a quartic Hermite interpolant p that satisfies:

p(i)=f _(i) , p*(i)=d _(i) and p**(i)=κ_(i)  (2)

[0044] where the components of p are polynomials on each parameter interval [i, i+1]. Also as above, the interpolant advantageously preserves the convexity of the data, that is, under the assumption:

sign(κ_(i))=sign(d _(i)×(f _(i+1) −f ₁))=sign(f _(i) −f _(i-1))×d _(i)), for all i  (3)

[0045] the curvature of the interpolant must have the same sign throughout.

[0046] According to one embodiment of the present invention, p is constructed on an interval-by-interval basis. Also, it should be noted that by constructing p on an interval-by-interval basis, the resulting interpolant p is necessarily twice continuously differentiable. The following will describe constructing the interpolant for the interval [0,1], but it should be understood that the interpolant can similarly be constructed for any of a number of other intervals without departing from the spirit and scope of the present invention. Turning attention to the interval [0,1], then, p can be written in its Bezier form as: $\begin{matrix} {{{p(t)} = {\sum\limits_{j = 0}^{4}\quad {b_{j}{B_{j}(t)}}}},} & \left( {0 \leq t \leq 1} \right) \end{matrix}$

[0047] where B_(j)(t) is defined as ${B_{j}(t)}:={\begin{pmatrix} 4 \\ j \end{pmatrix}{{t^{j}\left( {1 - t} \right)}^{({4 - j})}.}}$

[0048] As shown in FIG. 3, then, a method of modeling each of at least one section of a curve begins by providing positions (f₀, f₁) and associated directions (d₀, d₁) and curvatures (κ₀, κ₁), as shown in block 10. Then, points b₀, b₁, b₂, b₃ and b₄ can be identified, as shown in block 12. In this regard, points b₀, b₁, b₂, b₃ and b₄ can be identified based upon the provided positions (f₀, f₁) and associated directions (d₀, d₁) and curvatures (κ₀, κ₁). After points b₀, b₁, b₂, b₃ and b₄ have been identified, then, the quartic interpolant p(t) can be determined over an interval (0≦t≦1) based upon points b₀, b₁, b₂, b₃ and b₄ to thereby model the section of the curve, as illustrated in block 14. As will be appreciated, the interpolant p(t) has a direction and curvature equal to d_(i) and κ_(i), respectively, at t=i. Similarly, the interpolant p(t) has a direction and curvature equal to d_(i+1), and κ_(i+1), respectively, at t=i+1.

[0049] More particularly as to identifying points b₀, b₁, b₂, b₃ and b₄, the requirements of interpolation of position and tangent can result in relations among the coefficients b_(j), which can be written as follows:

b₀=f₀ and b₄=f₁

b ₁=(1−ρ₀)b ₀+ρ₀ {tilde over (b)} and b ₃=(1−ρ₁)b ₄+ρ₁ {tilde over (b)}  (5)

[0050] for some 0<ρ₀, ρ₁≦1, where ρ₀ and ρ₁ represent weights in writing b₁ and b₃ as convex combinations of b₀, b₄ and {tilde over (b)}. Further, as shown in equation (5), a control point {tilde over (b)} can be defined as a point at the intersection of a line through point f₀ parallel to direction d₀, and a line through point f₁ and parallel to direction d₁. More particularly, control point {tilde over (b)} can be defined as follows: $\overset{\sim}{b}:=\frac{{\left( {f_{1} \times d_{1}} \right)d_{0}} + {\left( {d_{0} \times f_{0}} \right)d_{1}}}{d_{0} \times d_{1}}$

[0051] According to embodiments of the present invention, b₂ is identified to preserve convexity in the interpolant. In this regard, b₂ is typically identified as a point within a region bounded by a triangular shape defined by points {tilde over (b)}, b₁ and b₃. As used herein, within a region bounded by a triangular shape, b₂ shall lie either within the interior or on the boundary of the triangular shape defined by points {tilde over (b)}, b₁ and b₃. More particularly, b₂ is typically identified such that the following condition is met to preserve convexity:

b ₂=α₀ b ₁+α₁ b ₃+(1−α₀+α₁){tilde over (b)}  (6)

[0052] for some α₀, α₁≧0 satisfying α₀+α₁≦1, as shown in FIG. 4. In this regard, α₀ and α₁ represent weights in writing b₂ as a convex combination of b₁, b₃ and {tilde over (b)}. Now, in order to meet the curvature requirements at the endpoints, κ_(i) may be determined as follows: $\begin{matrix} {\kappa_{0} = \frac{3d_{0} \times \left( {b_{2} - b_{1}} \right)}{4{{b_{1} - b_{0}}}^{2}}} & {and} & {\kappa_{1} = \frac{3\left( {b_{3} - b_{2}} \right) \times d_{1}}{4{{b_{4} - b_{3}}}^{2}}} \end{matrix}$

[0053] Substituting $\begin{matrix} \begin{matrix} {{b_{2} - b_{1}} = {{\left( {1 - \alpha_{0}} \right)\left( {\overset{\sim}{b} - b_{1}} \right)} + {\alpha_{1}\left( {b_{3} - \overset{\sim}{b}} \right)}}} \\ {\quad {= {{\left( {1 - \alpha_{0}} \right)\left( {1 - \rho_{0}} \right)\left( {\overset{\sim}{b} - b_{0}} \right)} + {{\alpha_{1}\left( {1 - \rho_{1}} \right)}\left( {b_{4} - \overset{\sim}{b}} \right)}}}} \end{matrix} \\ \begin{matrix} {{b_{3} - b_{2}} = {{\alpha_{0}\left( {\overset{\sim}{b} - b_{1}} \right)} + {\left( {1 - \alpha_{1}} \right)\left( {b_{3} - \overset{\sim}{b}} \right)}}} \\ {= {{{\alpha_{0}\left( {1 - \rho_{0}} \right)}\left( {\overset{\sim}{b} - b_{0}} \right)} + {\left( {1 - \alpha_{1}} \right)\left( {1 - \rho_{1}} \right)\left( {b_{4} - \overset{\sim}{b}} \right)}}} \end{matrix} \end{matrix}$

[0054] into the equations for κ₀ and κ₁ above, κ₀ and κ₁ can now be determined according to the following: $\begin{matrix} \begin{matrix} {\kappa_{0} = \frac{3{\alpha_{1}\left( {1 - \rho_{1}} \right)}d_{0} \times \left( {b_{4} - \overset{\sim}{b}} \right)}{4\rho_{0}^{2}{{\overset{\sim}{b} - b_{0}}}^{2}}} & {and} & {\kappa_{1} = \frac{3{\alpha_{0}\left( {1 - \rho_{0}} \right)}\left( {\overset{\sim}{b} - b_{0}} \right) \times d_{1}}{4\rho_{1}^{2}{{b_{4} - \overset{\sim}{b}}}^{2}}} \end{matrix} & (7) \end{matrix}$

[0055] Then, R₀ and R₁ can be defined as follows: $\begin{matrix} \begin{matrix} {R_{0}:=\frac{4\kappa_{0}{{\overset{\sim}{b} - b_{0}}}^{2}}{3d_{0} \times \left( {b_{4} - \overset{\sim}{b}} \right)}} & {and} & {R_{1}:=\frac{4\kappa_{1}{{b_{4} - \overset{\sim}{b}}}^{2}}{3\left( {\overset{\sim}{b} - b_{0}} \right) \times d_{1}}} \end{matrix} & (8) \end{matrix}$

[0056] With R₀ and R₁ defined according to equation (8), equation (7) can be solved for α₀ and α₁ such that α₀ and α₁ can be determined according to the following equation (9): $\begin{matrix} \begin{matrix} {\alpha_{0} = \frac{R_{1}\rho_{1}^{2}}{1 - \rho_{0}}} & {and} & {\alpha_{1} = \frac{R_{0}\rho_{0}^{2}}{1 - \rho_{1}}} \end{matrix} & (9) \end{matrix}$

[0057] Therefore, on each segment, f_(i) to f_(i+1), any value of ρ₀ and ρ₁ may be used, subject to the following conditions: $\begin{matrix} \begin{matrix} {{0 < \rho_{0}},{\rho_{1} < 1}} & {and} & {{\frac{R_{1}\rho_{1}^{2}}{1 - \rho_{0}} + \frac{R_{0}\rho_{0}^{2}}{1 - \rho_{1}}} \leq 1} \end{matrix} & (10) \end{matrix}$

[0058] Then, if the Bezier coefficients are prescribed as above, the resultant spline p will be C²-continuous (with respect to arclength) and will interpolate p(t) over the interval [0, 1], as in equation (2).

[0059] It should be noted that the conditions of (10) always have infinitely many solutions. For example, the sum in the second condition is equal to one when the following is true: $\begin{matrix} {\rho_{0} = {\rho_{1} = {\rho_{cnt} = \frac{\sqrt{1 + {4\left( {R_{0} + R_{1}} \right)}} - 1}{2\left( {R_{0} + R_{1}} \right)}}}} & (11) \end{matrix}$

[0060] And since 0<ρ_(cnt)<1, and the sum in the second condition of (10) is increasing in ρ₀ and ρ₁, any choice satisfying 0<ρ₀, ρ₁≦ρ_(cnt) can result in a valid interpolant p. It should be noted that these are not the only possible solutions and, as such, should not be taken to limit the scope of the present invention.

[0061] As described above, one aspect of the present invention provides a method of interpolating a position, tangent and curvature by a twice continuously differentiable quartic spline. In one advantageous embodiment, the values of ρ₀ and ρ₁ are defined such that the interpolant p is sixth-order accurate. According to this embodiment, suppose f comprises a smooth parametric curve with non-vanishing curvature, and t, is an increasing sequence of parameter values. Let f_(i), d_(i) and κ_(i) be given as in equation (1) above, and define the following: $\begin{matrix} {h_{i}:={{f_{i + 1} - f_{i}}}} & {and} & {h:={\max\limits_{i}\quad h_{i}}} \end{matrix}$

[0062] where h_(i) represents local step sizes or the distances between adjacent points, and h represents the maximum distance between any two adjacent points over all of i. For each segment, f_(i) to f_(i+1), then, ρ₀ and ρ₁ can be defined to achieve sixth-order accuracy of interpolant p(t) if the conditions in (10) are met, and the following equation (12) is satisfied: $\begin{matrix} {\rho_{0} = {\rho_{1} = {\frac{1}{2} + {A_{i}h_{i}^{2}} + {O\left( h_{i}^{3} \right)}}}} & (12) \end{matrix}$

[0063] where A_(i) is bounded uniformly in h. Also in equation (12), the term O, often referred to as “big-O,” represents a measure of the difference between ρ₀ or ρ₁ and the quantity ½+A_(i)h_(i) ², as such is defined by how much closer ρ₀ or ρ₁ gets to ½+A_(i)h_(i) ² as the number of known data points on f (i.e., f_(i), f_(i+1), etc.) increase. As such, the term O(h_(i) ³) in equation (12) indicates that ρ₀ and ρ₁ get g³ times closer to ½+A_(i)h_(i) ² if the number of known data points on f increases by a factor of g. By so defining ρ₀ and ρ₁, the following conclusion holds:

dist(f,p)=O(h ⁶).

[0064] In other terms, by defining p₀ and p₁ to satisfy (10) and (12), the distance between f and the interpolant p(t) is sixth order accurate such that if the number of known data points of f increase by a factor of g, the interpolant will get g⁶ times closer to f. The distance between f and the interpolant p(t) can be defined in any of a number of different known manners, including according to the Hausdorff technique.

[0065] To illustrate that defining ρ₀ and ρ₁ as in (12) results in sixth order accuracy of the interpolant p, assume that the curve f is smooth with |κ|≧c>0, where c represents a measure of curvature bounded away from zero. Parameterizing f by arclength s, then, results in the following: ${f^{\prime}(s)}:=\begin{bmatrix} {\cos \quad {\theta (s)}} \\ {\sin \quad {\theta (s)}} \end{bmatrix}$

[0066] with θ defined as the indefinite integral of the curvature,

θ′(s)=κ(s).

[0067] Again considering only one interval, designated f₀ to f₁, it can be assumed without loss of generality that f(0) is halfway along the arc from f₀ to f ₁ so that f₀=f(−ε) and f₁=f(ε) for some ε>0. Then, If |f₁−f₀|=2ε+I(ε³) that defining ρ₀ and ρ₁ to satisfy (10) and (12) results indist(f,p_(f))=O(ε⁶). It can also be assumed without loss of generality that θ(0)=0, i.e., that ${f^{''}(0)} = {\begin{bmatrix} 1 \\ 0 \end{bmatrix}.}$

[0068] Under such assumptions, the second coordinates of the curve f and the interpolant p can be considered as functions of their first coordinates by the implicit function theorem (for sufficiently small ε). That is, each curve can be parameterized by its first coordinate: $\begin{matrix} {p:\left. \xi\mapsto\begin{bmatrix} \xi \\ {y_{p}(\xi)} \end{bmatrix} \right.} & {and} & {f:\left. \xi\mapsto{\begin{bmatrix} \xi \\ {y_{f}(\xi)} \end{bmatrix}.} \right.} \end{matrix}$

[0069] From this point of view, the function y_(p) matches y_(f) to third order at two values of ξ that are O(ε) apart. As such, by the standard error estimate for two-point Hermite interpolation, the error is of order O(ε⁶) provided that d⁶y_(p)/dξ⁶ is bounded uniformly in ε.

[0070] Considering the Bezier form of the quartic polynomial p, ${{p(t)} = {\begin{bmatrix} {x(t)} \\ {y(t)} \end{bmatrix} = {\sum\limits_{j = 0}^{4}{b_{j}{B_{j}(t)}}}}},$

[0071] d^(n)y_(p)/dξ^(n) is a linear combination of the functions: $\left\{ {{\frac{y^{(i)}{\prod\limits_{\upsilon = 1}^{n - 1}x^{(j_{\upsilon})}}}{\left( x^{\prime} \right)^{{2n} - 1}}:{1 \leq i}},{{j_{\upsilon} \leq 4};{{i + {\sum\limits_{\upsilon = 1}^{n - 1}j_{\upsilon}}} = {{2n} - 1}}}} \right\}$

[0072] for any n≧1, as can be verified by inductively applying the chain rule, as such is known to those skilled in the art. In particular, the terms in d⁶y_(p)/dξ⁶ are of the form: $\frac{y^{(i)}{\prod\limits_{\upsilon = 1}^{5}x^{(j_{\upsilon})}}}{\left( x^{\prime} \right)^{11}}$

[0073] where i+Σ_(υ)j_(υ)=11 and 1≦i, j_(υ)≦4. The boundedness of these expressions can then be determined from the following assertions: $\begin{matrix} \begin{matrix} {{p^{\prime}(t)} = {\begin{bmatrix} {\frac{1}{2}ɛ} \\ 0 \end{bmatrix} + {O\left( ɛ^{2} \right)}}} & {and} & {{{{p^{(i)}(t)}} = {{{O\left( ɛ^{i} \right)}\quad {for}\quad i} = 2}},3,4} \end{matrix} & (13) \end{matrix}$

[0074] To illustrate how the assertions hold true, consider that the Bezier coefficients of the derivatives of p are as follows:

[0075] p′: 4(b₁−b₀) 4(b₂−b₁) 4(b₃−b₂) 4(b₄−b₃)

[0076] p″: 12(b₂−2b₁+b₀) 12(b₃−2b₂+b₃) 12(b₄−2b₃+b₂)

[0077] p′″: 24(b₃−3b₂+3b₁−b₀) 24(b₃−3b₂+3b₁−b₀)

[0078] p″: 24(b₄−4b₃+6b₂−4b₁+b₀)

[0079] Each of the Bezier coefficients can be expressed in terms of:

υ₀ {tilde over (b)}−b ₀ and υ₁ =b ₄ −{tilde over (b)}

[0080] by noting that:

b ₁ −b ₀=ρυ₀,

b ₂ −b ₁=(1−ρp)[(1−α₀)υ₀+α_(1υ) ₁],

b ₃ −b ₂=(1−ρ)[α₀υ₀+(1−α₁)υ₁] and

b ₄ −b ₃ =ρυ ₁.

[0081] Also, from the definition of {tilde over (b)}, it can be shown that: $\begin{matrix} {\upsilon_{0} = {\frac{a \times d_{1}}{d_{0} \times d_{1}}d_{0}}} & {and} & {{\upsilon_{1} = {\frac{d_{0} \times a}{d_{0} \times d_{1}}d_{1}}},} \end{matrix}$

[0082] where α:=f₁−f₀.

[0083] To expand the coefficients of p^((i)) in terms of ε, θ is initially expanded to the third order at s=0,

θ(s)=θ₁ s+θ ₂ s ²+θ₃ s ³ +O(s ⁴),

[0084] where θ_(k):=θ^((k))(0)/k! is independent of s. Then, from the definitions, it can be shown that: $\begin{matrix} {{d_{0} = {\begin{bmatrix} {\cos \quad {\theta \left( {- ɛ} \right)}} \\ {\sin \quad {\theta \left( {- ɛ} \right)}} \end{bmatrix} = {\begin{bmatrix} {1 - {\frac{\theta_{1}^{2}}{2}ɛ^{2}} + {\theta_{1}\theta_{2}ɛ^{3}}} \\ {{{- \theta_{1}}ɛ} + {\theta_{2}ɛ^{2}} - {\left( {\theta_{3} - {\theta_{1}^{3}/6}} \right)ɛ^{3}}} \end{bmatrix} + {O\left( ɛ^{4} \right)}}}},} \\ {{d_{1} = {\begin{bmatrix} {\cos \quad {\theta (ɛ)}} \\ {\sin \quad {\theta (ɛ)}} \end{bmatrix} = {\begin{bmatrix} {1 - {\frac{\theta_{1}^{2}}{2}ɛ^{2}} - {\theta_{1}\theta_{2}ɛ^{3}}} \\ {{\theta_{1}ɛ} + {\theta_{2}ɛ^{2}} + {\left( {\theta_{3} - {\theta_{1}^{3}/6}} \right)ɛ^{3}}} \end{bmatrix} + {O\left( ɛ^{4} \right)}}}},} \\ {a = {{\int_{- ɛ}^{ɛ}{\begin{bmatrix} {\cos \quad {\theta (\sigma)}} \\ {\sin \quad {\theta (\sigma)}} \end{bmatrix}{\sigma}}} = {\begin{bmatrix} {{2\quad ɛ} - {\frac{\theta_{1}^{2}}{3}ɛ^{3}}} \\ {\frac{2\quad \theta_{2}}{3}ɛ^{3}} \end{bmatrix} + {{O\left( ɛ^{5} \right)}.}}}} \end{matrix}$

[0085] Thereafter, the following expansions can be determined: $\begin{matrix} {{{d_{0} \times d_{1}} = {{2\quad \theta_{1}ɛ} + {\left( {{2\quad \theta_{3}} - \frac{4\quad \theta_{1}^{3}}{3}} \right)ɛ^{3}} + {O\left( ɛ^{5} \right)}}},} \\ {{d_{0} \times a} = {{2\quad \theta_{1}ɛ^{2}} - {\frac{4}{3}\theta_{2}ɛ^{3}} + {\left( {{2\quad \theta_{3}} - {\frac{2}{3}\theta_{1}^{3}}} \right)ɛ^{4}} + {{O\left( ɛ^{5} \right)}\quad {and}}}} \\ {{a \times d_{1}} = {{2\quad \theta_{1}ɛ^{2}} + {\frac{4}{3}\theta_{2}ɛ^{3}} + {\left( {{2\quad \theta_{3}} - {\frac{2}{3}\theta_{1}^{3}}} \right)ɛ^{4}} + {{O\left( ɛ^{5} \right)}.}}} \end{matrix}$

[0086] Next, the vectors υ₀ and υ₁ can be determined as follows: $\begin{matrix} {{\upsilon_{0} = {{\frac{a \times d_{1}}{d_{0} \times d_{1}}d_{0}} = {\begin{bmatrix} {ɛ + {\frac{2\quad \theta_{2}}{3\quad \theta_{1}}ɛ^{2}} - {\frac{\theta_{1}^{2}}{6}ɛ^{3}}} \\ {{{- \theta_{1}}ɛ^{2}} + {\frac{\theta_{2}}{3}ɛ^{3}}} \end{bmatrix} + {O\left( ɛ^{4} \right)}}}},} \\ {\upsilon_{1} = {{\frac{d_{0} \times a}{d_{0} \times d_{1}}d_{1}} = {\begin{bmatrix} {ɛ - {\frac{2\quad \theta_{2}}{3\quad \theta_{1}}ɛ^{2}} - {\frac{\theta_{1}^{2}}{6}ɛ^{3}}} \\ {{\theta_{1}ɛ^{2}} + {\frac{\theta_{2}}{3}ɛ^{3}}} \end{bmatrix} + {{O\left( ɛ^{4} \right)}.}}}} \end{matrix}$

[0087] It will be noted that in the above equations for υ₀ and υ₁, it is presumed that κ≠0, as otherwise the numerator vanishes for some value of s ∈[−ε,ε].

[0088] Now, since κ(S)=θ′(S), κ₀ and κ₁ can be determined as follows:

κ₀=κ(−ε)=θ₁−2θ₂ε+3θ₃ε³ +O(ε³), and

κ₁=κ(ε)=θ₁+2θ₂ε+3θ₃ε³ +O(ε³).

[0089] Based upon the definitions of R₀ and R₁ in equation (8) and the foregoing expressions for κ₀ and κ₁, R₀ and R₁ can be rewritten as follows: $\begin{matrix} {{R_{0} = {{\frac{4}{3}\frac{{\kappa_{0}\left( {a \times d_{1}} \right)}^{2}}{\left( {d_{0} \times a} \right)\left( {d_{0} \times d_{1}} \right)^{2}}} = {\frac{2}{3} + \frac{2\left( {{9\quad \theta_{1}^{4}} - {20\quad \theta_{2}^{2}} + {18\quad \theta_{1}\theta_{3}}} \right)}{27\quad \theta_{1}^{2}} + {O\left( ɛ^{3} \right)}}}},} \\ {R_{1} = {{\frac{4}{3}\frac{{\kappa_{1}\left( {d_{0} \times a} \right)}^{2}}{\left( {a \times d_{1}} \right)\left( {d_{0} \times d_{1}} \right)^{2}}} = {\frac{2}{3} + \frac{2\left( {{9\quad \theta_{1}^{4}} - {20\quad \theta_{2}^{2}} + {18\quad \theta_{1}\theta_{3}}} \right)}{27\quad \theta_{1}^{2}} + {{O\left( ɛ^{3} \right)}.}}}} \end{matrix}$

[0090] Combining the foregoing equations for R₀ and R₁ with the definition of ρ₀ and ρ₁ from equation (12) gives: $\alpha_{i} = {\frac{R_{i}\rho}{1 - \rho} = {\frac{1}{3} + {\left( {{2A_{0}} + \frac{\theta_{1}^{2}}{3} - \frac{20\quad \theta_{2}^{2}}{27} + \frac{\theta_{3}}{3\quad \theta_{1}}} \right)ɛ^{2}} + {O\left( ɛ^{3} \right)}}}$

[0091] In turn, the difference between consecutive points b_(i) can be determined as follows: $\begin{matrix} {{{b_{1} - b_{0}} = {{\rho \quad \upsilon_{0}} = {\begin{bmatrix} {{\frac{1}{2}ɛ} + {\frac{\theta_{2}}{3\quad \theta_{1}}ɛ^{2}} + {\left( {A_{0} - \frac{\theta_{1}^{2}}{12}} \right)ɛ^{3}}} \\ {{{- \frac{\theta_{1}}{2}}ɛ^{2}} + {\frac{\theta_{2}}{6}ɛ^{3}}} \end{bmatrix} + {O\left( ɛ^{4} \right)}}}},} \\ {{b_{2} - b_{1}} = {\left( {1 - \rho} \right)\left\lbrack {{\left( {1 - \alpha_{0}} \right)\upsilon_{0}} + {\alpha_{1}\upsilon_{1}}} \right\rbrack}} \\ {\quad {{= {\begin{bmatrix} {{\frac{1}{2}ɛ} + {\frac{\theta_{2}}{9\quad \theta_{1}}ɛ^{2}} - {\left( {A_{0} + \frac{\theta_{1}^{2}}{12}} \right)ɛ^{3}}} \\ {{{- \frac{\theta_{1}}{6}}ɛ^{2}} + {\frac{\theta_{2}}{6}ɛ^{3}}} \end{bmatrix} + {O\left( ɛ^{4} \right)}}},}} \\ {{b_{3} - b_{2}} = {\left( {1 - \rho} \right)\left\lbrack {{\alpha_{0}\upsilon_{0}} + {\left( {1 - \alpha_{0}} \right)\upsilon_{1}}} \right\rbrack}} \\ {\quad {= {\begin{bmatrix} {{\frac{1}{2}ɛ} - {\frac{\theta_{2}}{9\quad \theta_{1}}ɛ^{2}} - {\left( {A_{0} + \frac{\theta_{1}^{2}}{12}} \right)ɛ^{3}}} \\ {{\frac{\theta_{1}}{6}ɛ^{2}} + {\frac{\theta_{2}}{6}ɛ^{3}}} \end{bmatrix} + {{O\left( ɛ^{4} \right)}\quad {and}}}}} \\ {{b_{4} - b_{3}} = {{\rho \quad \upsilon_{1}} = {\begin{bmatrix} {{\frac{1}{2}ɛ} - {\frac{\theta_{2}}{3\quad \theta_{1}}ɛ^{2}} + {\left( {A_{0} - \frac{\theta_{1}^{2}}{12}} \right)ɛ^{3}}} \\ {{\frac{\theta_{1}}{2}ɛ^{2}} + {\frac{\theta_{2}}{6}ɛ^{3}}} \end{bmatrix} + {{O\left( ɛ^{4} \right)}.}}}} \end{matrix}$

[0092] From the foregoing equations, it follows that ${{b_{j + 1} - b_{j}} = {\begin{bmatrix} {\frac{1}{2}ɛ} \\ 0 \end{bmatrix} + {O\left( ɛ^{2} \right)}}},$

[0093] |b_(j+2)−2b_(j+1)+b_(j)|=I(ε²),

|b _(j+1)−3b _(j+2)+3b _(j+1) −b _(j) |=O(ε³) and |b _(j+4)−4b _(j+3)+6b _(j+2)−4b _(j+1) +b _(j) |=O(ε⁴).

[0094] As shown from the foregoing, then, the foregoing equations satisfy the assertions in (13) above.

[0095] Now suppose that the data came from endpoints of an arc of angle θ<π of some circle of curvature κ. According to another embodiment of the present invention, to preserve the symmetry of such a situation, ρ₀ and ρ₁ can be defined such that ρ₀=ρ₁=:ρ. The heuristic choice of a favorable ρ can force the curvature of the interpolant to be κ at the midpoint, that is, p**(½)=κ_(i). In terms of the spline coefficients, one can verify that $\begin{matrix} {{p^{\prime}\left( {1/2} \right)} = \frac{b_{4} + {2b_{3}} - {2b_{1}} - b_{0}}{2}} & {and} & {{p^{''}\left( {1/2} \right)} = {3{\left( {b_{4} - {2b_{2}} + b_{0}} \right).}}} \end{matrix}$

[0096] For any fixed choice of ρ₀ and ρ₁, the construction devised above is scale-invariant (provided the curvature is scaled by 1/s to correspond to a scaling by s in the relative positions). So without loss of generality, it can be assumed that κ=1. The construction is also invariant under translation and rotation, so it can be assumed that

[0097] f₀=(0,0), d₀=(1,0), f₁=(sinθ,1−cosθ), d₁=(cosθ,sinθ) and κ₀=κ₁=1.

[0098] In such a case, as shown in FIG. 5, b₀=(0,0), b₄=(sinθ,1−cosθ) and {tilde over (b)}=(sinθ/(1+cosθ),0).

[0099] From equations (8) and (9) above, R₀=R, =4/3(1+cosθ) and $\alpha_{0} = {\alpha_{1} = {\alpha:={\frac{4\quad \rho^{2}}{3\left( {1 - \rho} \right)\left( {1 + {\cos \quad \vartheta}} \right)}.}}}$

[0100] Defining υ₀:={tilde over (b)}−b₀ and υ₁:={tilde over (b)}−b₄, b₁ and b₃ can be determined as follows: b₁=b₀+ρυ₀, and b₃=b₄+ρυ₁. And since υ₀−υ₁=b₄−b₀, b₄+2b₃−2b₁−b₀=(3−2ρ)(υ₀−υ₁).

[0101] From equation (6), then, b₂=b−α(1−ρ)(υ₀+υ₁). Therefore, it can be shown that: ${b_{4} - {2b_{2}} + b_{0}} = {{{- \left( {1 - {2\quad {\alpha \left( {1 - \rho} \right)}}} \right)}\left( {\upsilon_{0} + \upsilon_{1}} \right)} = {{- \left( {1 - \frac{8\quad \rho^{2}}{3\left( {1 + {\cos \quad \vartheta}} \right)}} \right)}{\left( {\upsilon_{0} + \upsilon_{1}} \right).}}}$

[0102] Consequently, the first and second derivatives for the interpolant p at t=½ can be determined as follows: $\begin{matrix} {{p^{\prime}\left( {1/2} \right)} = \frac{\left( {3 - {2\quad \rho}} \right)\left( {\upsilon_{0} - \upsilon_{1}} \right)}{2}} & {and} & {{p^{''}\left( {1/2} \right)} = {{- \left( {3 - \frac{8\quad \rho^{2}}{1 + {\cos \quad \vartheta}}} \right)}{\left( {\upsilon_{0} + \upsilon_{1}} \right).}}} \end{matrix}$

[0103] Since |υ₀|=|υ₁|, it can also be shown that υ₀−υ₁ and υ₀+υ₁ are orthogonal, Therefore, (υ₀−υ₁)×(υ₀+υ₁)=−|υ₀−υ₁|·|υ₀+υ₁|. It can further be shown that $\begin{matrix} {{p^{**}\left( {1/2} \right)} = \frac{{p^{\prime}\left( {1/2} \right)} \times {p^{''}\left( {1/2} \right)}}{{{p^{\prime}\left( {1/2} \right)}}^{3}}} \\ {= \frac{\left( \frac{3 - {2\quad \rho}}{2} \right)\left( {3 - \frac{8\quad \rho^{2}}{1 + {\cos \quad \vartheta}}} \right){{\upsilon_{0} - \upsilon_{1}}}{{\upsilon_{0} + \upsilon_{1}}}}{\frac{1}{8}\left( {3 - {2\quad \rho}} \right)^{3}{{\upsilon_{0} - \upsilon_{1}}}^{3}}} \\ {= {\frac{4\left( {3 - \frac{8\quad \rho^{2}}{1 + {\cos \quad \vartheta}}} \right){{\upsilon_{0} + \upsilon_{1}}}}{\left( {3 - {2\quad \rho}} \right)^{2}{{\upsilon_{0} - \upsilon_{1}}}^{2}}.}} \end{matrix}$

[0104] By noting that $\rho \neq \frac{3}{2}$

[0105] and substituting in $\begin{matrix} {{{\upsilon_{0} - \upsilon_{1}}} = \sqrt{2\left( {1 - {\cos \quad \vartheta}} \right)}} & \begin{matrix} {and} & {{{{\upsilon_{0} + \upsilon_{1}}} = \frac{2\left( {1 - {\cos \quad \vartheta}} \right)}{\sqrt{2\left( {1 + {\cos \quad \vartheta}} \right)}}},} \end{matrix} \end{matrix}$

[0106] the above equation for p**(½) can be simplified to the following quadratic equation when p**(½) is set equal to 1: ${{\left( {\sqrt{2\left( {1 - {\cos \quad \vartheta}} \right)} + \frac{8}{1 + {\cos \quad \vartheta}}} \right)\rho^{2}} - {3\sqrt{2\left( {1 + {\cos \quad \vartheta}} \right)}\rho} + \left( {{\frac{9}{4}\sqrt{2\left( {1 + {\cos \quad \vartheta}} \right)}} - 3} \right)} = 0.$

[0107] The larger solution of the foregoing quadratic equation can then be denoted by ρ_(opt) (θ), where $0 < {\rho_{opt}(\vartheta)} \leq \frac{1}{2}$

[0108] for all −π<θ<π.

[0109] It is advantageous that ρ₀ ρ₁=ρ_(opt) (θ) for circular data as above, and that the resulting interpolant provides sixth order accuracy in general. Therefore, ρ₀ and ρ₁ can be defined to continuously depend on the data and meet the circular data requirement (as defined in equation 12). Therefore, for general data, it can be shown that $\begin{matrix} {\rho_{0} = {\rho_{1} = {\rho \text{:}{= \rho_{crit} \cdot \left( \frac{\rho_{opt}(\vartheta)}{\rho_{crit}(\vartheta)} \right)}}}} & (14) \end{matrix}$

[0110] where θ is the angle between d₀ and d₁, and ρ_(cnt) (θ) is defined as ρ_(cnt) for a circular arc of angle θ. More particularly, ${\rho_{crit}(\vartheta)} = \frac{3\left( \sqrt{1 + \frac{32}{3\left( {1 + {\cos \quad \vartheta}} \right)} - 1} \right)\left( {1 + {\cos \quad \vartheta}} \right)}{16}$

[0111] (see equation (11)). The preceding equation results in a favorable choice of ρ for circular data and necessarily satisfies 0<ρ₀, ρ₁<ρ_(cnt) (see FIG. 6 for a plot of ρ_(opt) (θ)/ρ_(cnt) (θ) for −π<θ<π).

[0112] As indicated above, embodiments of the present invention provide an interpolant with sixth order accuracy. As will be appreciated, many different techniques for choosing ρ can provide the same accuracy, including the particularly simple choice ρ₀=ρ₁=min{{fraction (1/2)}, ρ_(cnt)}. Empirically, such a choice tends to produce overly flat curves for sparse data sets. The fact that the interpolant is sixth-order accurate follows from the above, and the claim that, supposing f is a smooth parameteric curve with non-vanishing curvature, then ρ as defined in equation (14) satisfies the following equation (15): $\begin{matrix} {\rho = {\frac{1}{2} + {Ah}^{2} + {O\left( h^{3} \right)}}} & (15) \end{matrix}$

[0113] for some A ∈1R independent of h.

[0114] To illustrate how the interpolant is sixth-order accurate by demonstrating that ρ satisfies equation (15), consider that cos θ=<d₀,d₁>, the inner product of the vectors d₀ and d₁. As such, from the expansions d₀×d₁, d₀×α and α×d₁ described above,

cosθ=1−2θ ₁ ² h ² +O(h ⁴).

[0115] Now, letting a, b and c denote the coefficients in the quadratic equation that defines ρ_(opt) (θ), it can be shown that:

[0116] a=6+3θ₁ ²h²+O(h⁴), b=−6+30₁ ²h²+O(h⁴), and $c = {\frac{3}{2} - {\frac{9}{4}\theta_{1}^{2}h^{2}} + {{O\left( h^{4} \right)}.}}$

[0117] Therefore, ${\rho_{opt}(\vartheta)} = {\frac{\sqrt{b^{2} - {4{ac}}} - b}{2a} = {\frac{1}{2} - {\frac{1}{2}\theta_{1}^{2}h^{2}} + {{O\left( h^{4} \right)}.}}}$

[0118] As above, $R_{i} = {\frac{2}{3} + {\frac{2\left( {{9\quad \theta_{1}^{4}} - {20\quad \theta_{2}^{2}} + {18\quad \theta_{1}\theta_{3}}} \right)}{27\quad \theta_{1}^{2}}h^{2}} + {O\left( h^{4} \right)}}$

[0119] for any curve, so ρ_(cnt) can be defined as: $\rho_{crit} = {\frac{\sqrt{57} - 3}{8} + {{O\left( h^{2} \right)}.}}$

[0120] Since θ=O(h) for any fixed sufficiently smooth curve,

ρ_(cnt)/ρ_(cnt)(θ)=1+O(h²).

[0121] Therefore, continuing with a more rigorous examination, it can be shown that $\rho_{i} = {\frac{1}{2} + {{O\left( h^{2} \right)}.}}$

[0122] As described above, points b₀, b₁, b₂, b₃ and b₄ are based upon the pair of positions (f₀, f₁) and associated directions (d₀, d₁) and curvatures (κ₀, κ₁). In turn, the interpolant p(t) is determined over the interval (0≦t≦1) based upon points b₀, b₁, b₂, b₃ and b₄ to thereby model a respective section of the curve. As such, it will be appreciated that the interpolant depends indirectly on the positions (f₀, f₁) and associated directions (d₀, d₁) and curvatures (κ₀, κ₁). Advantageously, then, a change in the interpolant p(t) is directly proportional to a change in the positions (f_(i), f_(i+1)), directions (d_(i), d_(i+1)) and curvatures (κ_(i), κ_(i+1)). In other terms, the interpolant depends continuously on the data (i.e., the positions, directions and curvatures). As such, small changes in any of the data (positions, directions and/or curvatures) will not produce large changes to the interpolant.

[0123] As will also be appreciated, the system, method and computer program product of embodiments of the present invention can have any of a number of different applications. For example, embodiments of the present invention can be utilized to provide a parametric family of curves for use in the design of aircraft, rotorcraft, spacecraft, automobiles, maritime vehicles, or the like. Thus, according to one typical scenario, the system, method and computer program product of embodiments of the pesent invention can be utilized to design an airfoil along a wing of a prospective aircraft. Initially, then, the airfoil can be specified by a plurality of data points, includeing a position, direction and curvature for each point. Thus, presuming the airfoil is specified by ten points, the points according to one exemplar embodiment may appear as in FIG. 7A and Table 1, below, with the curve fit according to embodiments of the present invention shown in FIG. 7B. TABLE 1 i f_(i) d₁ k₁ 0 (1.00, 0.00322) (−0.756, 0.655) 0.0677 1 (0.552, 0.337) (−0.8577, 0.5141) 0.311 2 (0.220, 0.468) (−1.00, 0.00) 4.32 3 (0.0617, 0.350) (−0.433, −0.901) 3.16 4 (0.00, 0.00) (0.00, −1.00) 0.391 5 (0.0676, −0.403) (0.537, −0.844) 4.73 6 (0.316, −0.539) (1.00, 0.00) 5.43 7 (0.514, −0.451) (0.796, 0.605) 1.19 8 (0.728, −0.243) (0.676, 0.737) 0.00 9 (1.00, 0.00) (0.942, 0.336) −5.23

[0124] Each data point can be represented by a circular arc of the appropriate curvature centered at the appropriate position and tangent to the appropriate direction vector at the respective point. The curve defined by the data points and the interpolation technique of embodiments of the present invention can be considered to be a member of a 40-parameter family of curves. In this regard, the parameters consist of the position (x and y coordinates), direction (a single parameter since |d_(i)|=1), and curvature at each point. In fact, not all parameters are active in such a case. Namely, the following parameters can be held in position: x₀=x₉=1, x₄=y₀=y₄=y₉=0, d₂=(−1,0), d₄=(0, −1), d₆=(1, 0) and κ₈=0. As such, the number of parameters can be reduced from 40 to 30. It will be noted that it has been found that such a 30-parameter family of curves is sufficiently robust to model airfoils for a wide variety of aircraft. While advantageous to accurately model airfoils, the system, method and computer program product can be utilized to model a wide variety of other surfaces for other applications, including automotive applications and the like.

[0125] In another typical scenario, the curve $\begin{matrix} {{f(t)} = \begin{bmatrix} t \\ e^{t} \end{bmatrix}} & \left( {0 \leq t \leq 1} \right) \end{matrix}$

[0126] is interpolated to empirically demonstrate the sixth order accuracy of the interpolation technique. The results of such interpolation are shown in Table 2, below. TABLE 2 Points ρ_(min) ρ_(max) Error e_(i) $\frac{e_{i - 1}}{e_{i}}$

$\log_{2}\frac{e_{i - 1}}{e_{i}}$

2 .508499 .508499 9.5738 × 10⁻⁶ 3 .502127 .502150 1.6097 × 10⁻⁷ 59.476 5.8942 5 .500531 .500539 2.5669 × 10⁻⁹ 62.710 5.9706 9 .500133 .500135 4.0367 × 10⁻¹¹ 63.589 5.9907

[0127] As shown in Table 2, the first column specifies the number of knots interpolated at, where the knots were distributed uniformly in t. The next two columns specify the minimum and maximum ρ used on any of the intervals. The fourth column is the maximum error over the interval. Finally, the fifth and sixth columns display (or at least suggest) the sixth-order rate of convergence. Note that $\rho = {\frac{1}{2} + {O\left( h^{2} \right)}}$

[0128] was shown above.

[0129] Therefore, embodiments of the present invention are capable of modeling each section of a curve f to exactly match the positions, directions and curvatures for data provided for the respective section. Embodiments of the present invention are also capable of modeling each section of the curve to thereby preserve the shape properties of the curve f, such as by maintaining the convexity in the modeled section. In addition, by modeling sections of a curve, embodiments of the present invention can advantageously model the curve such that changes in the provided data only affect the curve locally. Embodiments of the present invention are further capable of modeling each section such that the model depends continuously on the data, and such that the model has sixth-order accuracy.

[0130] As shown in FIG. 8, the system of the present invention is typically embodied by a processing element and an associated memory device, both of which are commonly comprised by a computer 20 or the like. In this regard, as indicated above, the method of embodiments of the present invention can be performed by the processing element manipulating data stored by the memory device with any one of a number of commercially available computer software programs. The computer can include a display 22 for presenting information relative to performing embodiments of the method of the present invention, including the various distributions, models and/or conclusions as determined according to embodiments of the present invention. To plot information relative to performing embodiments of the method of the present invention, the computer can further include a printer 24.

[0131] Also, the computer 20 can include a means for locally or remotely transferring the information relative to performing embodiments of the method of the present invention. For example, the computer can include a facsimile machine 26 for transmitting information to other facsimile machines, computers or the like. Additionally, or alternatively, the computer can include a modem 28 to transfer information to other computers or the like. Further, the computer can include an interface (not shown) to a network, such as a local area network (LAN), and/or a wide area network (WAN). For example, the computer can include an Ethernet Personal Computer Memory Card International Association (PCMCIA) card configured to transmit and receive information to and from a LAN, WAN or the like.

[0132] According to one aspect of the present invention and as mentioned above, the system of the present invention, such as a processing element typically embodied by a computer, generally operates under control of a computer program product. The computer program product for performing the methods of embodiments of the present invention includes a computer-readable storage medium, such as the non-volatile storage medium, and computer-readable program code portions, such as a series of computer instructions, embodied in the computer-readable storage medium.

[0133] In this regard, FIG. 1 is a flowchart of methods, systems and program products according to the invention. It will be understood that each block or step of the flowchart, and combinations of blocks in the flowchart, can be implemented by computer program instructions. These computer program instructions may be loaded onto a computer or other programmable apparatus to produce a machine, such that the instructions which execute on the computer or other programmable apparatus create means for implementing the functions specified in the flowchart block(s) or step(s). These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart block(s) or step(s). The computer program instructions may also be loaded onto a computer or other programmable apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart block(s) or step(s).

[0134] Accordingly, blocks or steps of the flowchart support combinations of means for performing the specified functions, combinations of steps for performing the specified functions and program instruction means for performing the specified functions. It will also be understood that each block or step of the flowchart, and combinations of blocks or steps in the flowchart, can be implemented by special purpose hardware-based computer systems which perform the specified functions or steps, or combinations of special purpose hardware and computer instructions.

[0135] Many modifications and other embodiments of the invention will come to mind to one skilled in the art to which this invention pertains having the benefit of the teachings presented in the foregoing descriptions and the associated drawings. Therefore, it is to be understood that the invention is not to be limited to the specific embodiments disclosed and that modifications and other embodiments are intended to be included within the scope of the appended claims. Although specific terms are employed herein, they are used in a generic and descriptive sense only and not for purposes of limitation. 

What is claimed is:
 1. A method of modeling at least one section of a curve f, wherein modeling a section comprises: providing a pair of positions (f_(i), f_(i+1)) of the section of the curve, wherein the pair of positions includes associated directions (d_(i), d_(i+1)) and associated curvatures (κ_(i), κ_(i+1)); identifying points b₀, b₁, b₂, b₃ and b₄ based upon the pair of positions (f_(i), f_(i+1)) and associated directions (d_(i), d_(i+1)) and curvatures (κ_(i), κ_(i+1)); and determining a quartic interpolant p(t) over an interval (i≦t≦i+1) based upon points b₀, b₁, b₂, b₃ and b₄ to thereby model the section of the curve, wherein the interpolant p(t) has a position, direction and curvature equal to f_(i), d_(i) and κ_(i), respectively, at t=i, and wherein the interpolant p(t) has a position, direction and curvature equal to f_(i+1), d_(i+1), and κ_(i+1), respectively, at t=i+1.
 2. A method according to claim 1, wherein identifying points b₀, b₁, b₂, b₃ and b₄ comprises: defining a control point {tilde over (b)} based upon the pair of positions (f_(i), f_(i+1)) and the associated directions (d_(i), d_(i+1)); identifying a point b₁ and a point b₃, wherein point b₁ is identified on a segment {overscore (b₀{tilde over (b)})} and point b₃ is identified on a segment {overscore (b₄{tilde over (b)})}, and wherein b₀ equals f_(i) and b₄ equals f_(i+1); and identifying a point b₂ disposed within a region bounded by a triangular shape defined by points {tilde over (b)}, b₁ and b₂, wherein point b₂ is identified based upon the curvatures (κ_(i), κ_(i+1)).
 3. A method according to claim 1, wherein defining a control point {tilde over (b)} comprises defining a control point {tilde over (b)} as a point at the intersection of a line through point f₀ parallel to direction d₀, and a line through point f, and parallel to direction d₁.
 4. A method according to claim 3, wherein defining a control point comprises determining a control point according to the following: $\overset{\sim}{b}:={\frac{{\left( {f_{1} \times d_{1}} \right)d_{0}} + {\left( {d_{0} \times f_{0}} \right)d_{1}}}{d_{0} \times d_{1}}.}$


5. A method according to claim 2, wherein identifying a point b₁ and a point b₃ comprises: defining weights (ρ₀, ρ₁), wherein the weights are identified such that 0<ρ₀, ρ₁≦1; and identifying a point b₁ and a point b₃ according to the following: b ₁=(1−ρ₀)b ₀+ρ₀ {tilde over (b)} and b ₃=(1−ρ₁)b ₄+ρ₁ {tilde over (b)}.
 6. A method according to claim 2, wherein identifying a point b₂ comprises: defining weights (α₀, α₁) based upon the curvatures (κ_(i); κ_(i+1)), wherein the weights (α₀, α₁) are defined such that α₀, α₁≧0 and α₀+α₁≦1; and identifying point b₂ according to the following: b ₂ =α ₀ b ₁ +α ₁ b ₃+(1−α₀−α₁){tilde over (b)}.
 7. A method according to claim 6, wherein defining weights (α₀, α₁) comprises: defining weights (ρ₀, ρ₁) such that 0<ρ₀, ρ₁≦1; and defining weights (α₀, α₁) according to the following: $\begin{matrix} {\alpha_{0} = \frac{R_{1}\rho_{1}^{2}}{1 - \rho_{0}}} & {and} & {{\alpha_{1} = \frac{R_{0}\rho_{0}^{2}}{1 - \rho_{1}}},} \end{matrix}$

wherein $\begin{matrix} {R_{0}:=\frac{4\quad \kappa_{i}{{\overset{\sim}{b} - b_{0}}}^{2}}{3d_{i} \times \left( {b_{4} - \overset{\sim}{b}} \right)}} & {and} & {R_{1}:={\frac{4\quad \kappa_{i + 1}{{b_{4} - \overset{\sim}{b}}}^{2}}{3\left( {\overset{\sim}{b} - b_{0}} \right) \times d_{i + 1}}.}} \end{matrix}$


8. A method according to claim 1, wherein determining an interpolant p(t) comprises determining an interpolant p(t) according to the following: ${{p(t)} = {\sum\limits_{j = 0}^{4}{b_{j}{B_{j}(t)}}}},$

and wherein ${B_{j}(t)}:={\begin{pmatrix} 4 \\ j \end{pmatrix}{{t^{j}\left( {1 - t} \right)}^{({4 - j})}.}}$


9. A method according to claim 1, wherein identifying points b₀, b₁, b₂, b₃ and b₄ comprises identifying points b₀, b₁, b₂, b₃ and b₄ such that the quartic interpolant p(t) is sixth-order accurate over the interval (i≦t≦i+1).
 10. A method according to claim 1, wherein determining the quartic interpolant p(t) comprises determining the quartic interpolant p(t) such that a change in the interpolant p(t) is directly proportional to a change in at least one of the provided positions (f_(i), f_(i+1)), directions (d_(i), d_(i+1)) and curvatures (κ_(i), κ_(i+1)).
 11. A system for modeling at least one section of a curve f comprising: a processing element capable of receiving a pair of positions (f_(i), f_(i+1)) of a section of the curve to be modeled, wherein the pair of positions includes associated directions (d_(i), d_(i+1)) and associated curvatures (κ_(i), κ_(i+1)), wherein for each section the processing element is capable of identifying points b₀, b₁, b₂, b₃ and b₄ based upon the pair of positions (f_(i), f_(i+1)) and associated directions (d_(i), d_(i+1)) and curvatures (κ_(i), κ_(i+1)), wherein the processing element is also capable of determining a quartic interpolant p(t) over an interval (i≦t≦i+1) based upon points b₀, b₁, b₂, b₃ and b₄ to thereby model the respective section of the curve, wherein the interpolant p(t) has a position, direction and curvature equal to f_(i), d_(i) and κ_(i), respectively, at t=i, and wherein the interpolant p(t) has a position, direction and curvature equal to f_(i+1), d_(i+1), and κ_(i+1), respectively, at t=i+1.
 12. A system according to claim 11, wherein the processing element is capable of identifying points b₀, b₁, b₂, b₃ and b₄ by: defining a control point b based upon the pair of positions (f_(i), f_(i+1)) and the associated directions (d_(i), d_(i+1)); identifying a point b₁ and a point b₃, wherein point b₁ is identified on a segment {overscore (b₀{tilde over (b)})} and point b₃ is identified on a segment {overscore (b₄{tilde over (b)})}, and wherein b₀ equals f_(i) and b₄ equals f_(i+1); and identifying a point b₂ disposed within a region bounded by a triangular shape defined by points {tilde over (b)}, b₁ and b₂, wherein point b₂ is identified based upon the curvatures (κ_(i), κ_(i+1)).
 13. A system according to claim 12, wherein the processing element is capable of defining the control point {tilde over (b)} as a point at the intersection of a line through point f₀ parallel to direction d₀, and a line through point f₁ and parallel to direction d₁.
 14. A system according to claim 13, wherein the processing element is capable of defining the control point according to the following: $\overset{\sim}{b}:={\frac{{\left( {f_{1} \times d_{1}} \right)d_{0}} + {\left( {d_{0} \times f_{0}} \right)d_{1}}}{d_{0} \times d_{1}}.}$


15. A system according to claim 12, wherein the processing element is capable of identifying points b₁ and b₃ by: defining weights (ρ₀, ρ₁), wherein the weights are identified such that 0<ρ₀, ρ₁≦1; and identifying a point b₁ and a point b₃ according to the following: b ₁=(1−ρ₀)b ₀+ρ₀ {tilde over (b)} and b ₃=(1−ρ₁)b ₄+ρ₁ {tilde over (b)}.
 16. A system according to claim 12, wherein the processing element is capable of identifying point b₂ by: defining weights (α₀, α₁) based upon the curvatures (κ_(i), κ_(i+1)), wherein the weights (α₀, α₁) are defined such that α₀, α₁≧0 and α₀+α₁≦1; and identifying point b₂ according to the following: b ₂=α₀ b ₁+α₁ b ₃+(1−α₀−α₁){tilde over (b)}.
 17. A system according to claim 16, wherein the processing element is capable of defining weights (α₀, α₁) by: defining weights (ρ₀, ρ₁) such that 0<ρ₀, ρ₁≦1; and defining weights (α₀, α₁) according to the following: $\begin{matrix} {\alpha_{0} = \frac{R_{1}\rho_{1}^{2}}{1 - \rho_{0}}} & {and} & {{\alpha_{1} = \frac{R_{0}\rho_{0}^{2}}{1 - \rho_{1}}},} \end{matrix}$

wherein $\begin{matrix} {R_{0}:=\frac{4\kappa_{i}{{\overset{\sim}{b} - b_{0}}}^{2}}{3d_{i} \times \left( {b_{4} - \overset{\sim}{b}} \right)}} & {and} & {R_{1}:={\frac{4\kappa_{i + 1}{{b_{4} - \overset{\sim}{b}}}^{2}}{3\left( {\overset{\sim}{b} - b_{0}} \right) \times d_{i + 1}}.}} \end{matrix}$


18. A system according to claim 11, wherein the processing element is capable of determining the interpolant p(t) according to the following: ${{p(t)} = {\sum\limits_{j = 0}^{4}\quad {b_{j}{B_{j}(t)}}}},$

and wherein ${B_{j}(t)}:={\begin{pmatrix} 4 \\ j \end{pmatrix}{{t^{j}\left( {1 - t} \right)}^{({4 - j})}.}}$


19. A system according to claim 11, wherein the processing element is capable of identifying points b₀, b₁, b₂, b₃ and b₄ such that the quartic interpolant p(t) is sixth-order accurate over the interval (i≦t≦i+1).
 20. A system according to claim 11, wherein the processing element is capable of determining the quartic interpolant p(t) such that a change in the interpolant p(t) is directly proportional to a change in at least one of the provided positions (f_(i), f_(i+1)), directions (d_(i), d_(i+1)) and curvatures (κ_(i), κ_(i+1)).
 21. A computer program product for modeling at least one section of a curve f said computer program product comprising a computer-readable storage medium having computer-readable program code portions stored therein, the computer-readable program portions comprising: a first executable portion for receiving a pair of positions (f_(i), f_(i+1)) of a section of the curve to be modeled, wherein the pair of positions includes associated directions (d_(i), d_(i+1)) and associated curvatures (κ_(i), κ_(i+1)); a second executable portion for identifying points b₀, b₁, b₂, b₃ and b₄ based upon the pair of positions (f_(i), f_(i+1)) and associated directions (d_(i), d_(i+1)) and curvatures (κ_(i), κ_(i+1)); and a third executable portion for determining a quartic interpolant p(t) over an interval (i≦t≦i+1) based upon points b₀, b₁, b₂, b₃ and b₄ to thereby model the section of the curve, wherein the interpolant p(t) has a position, direction and curvature equal to f_(i), d_(i) and κ_(i), respectively, at t=i, and wherein the interpolant p(t) has a position, direction and curvature equal to f_(i+1), d_(i+1), and κ_(i+1), respectively, at t=i+1.
 22. A computer program product according to claim 21, wherein the second executable portion identifies points b₀, b₁, b₂, b₃ and b₄ by: defining a control point b based upon the pair of positions (f_(i), f_(i+1)) and the associated directions (d_(i), d_(i+1)); identifying a point b₁ and a point b₃, wherein point b₁ is identified on a segment {overscore (b₀{tilde over (b)})} and point b₃ is identified on a segment {overscore (b₄{tilde over (b)})}, and wherein b₀ equals f_(i) and b₄ equals f_(i+1); and identifying a point b₂ disposed within a region bounded by a triangular shape defined by points {tilde over (b)}, b₁ and b₂, wherein point b₂ is identified based upon the curvatures (κ_(i), κ_(i+1)).
 23. A computer program product according to claim 21, wherein the second executable portion defines the control point {tilde over (b)} as a point at the intersection of a line through point f₀ parallel to direction d₀, and a line through point f₁ and parallel to direction d₁.
 24. A computer program product according to claim 23, wherein the second executable portion defines the control point according to the following: $\overset{\sim}{b}:={\frac{{\left( {f_{1} \times d_{1}} \right)d_{0}} + {\left( {d_{0} \times f_{0}} \right)d_{1}}}{d_{0} \times d_{1}}.}$


25. A computer program product according to claim 22, wherein the second executable portion identifies points b₁ and b₃ by: defining weights (ρ₀, ρ₁) such that 0<ρ₀, ρ₁≦1; and identifying a point b₁ and a point b₃ according to the following: b ₁=(1−ρ₀)b ₀+ρ₀ {tilde over (b)} and b ₃=(1−ρ₁)b ₄+ρ₁ {tilde over (b)}.
 26. A computer program product according to claim 22, wherein the second executable portion identifies point b₂ by: defining weights (α₀, α₁) based upon the curvatures (κ_(i), κ_(i+1)), wherein the weights (α₀, α₁) are defined such that α₀, α₁≧0 and α₁+≦1; and identifying point b₂ according to the following: b ₂=α₀ b ₁+α₁ b ₃+(1−α₀−α₁){tilde over (b)}.
 27. A computer program product according to claim 26, wherein the second executable portion defines weights (α₀, α₁) by: defining weights (ρ₀, ρ₁), wherein the weights (ρ₀, ρ₁) are identified such that 0<ρ₀, ρ₁≦1; and defining weights (α₀, α₁) according to the following: $\begin{matrix} {\alpha_{0} = \frac{R_{1}\rho_{1}^{2}}{1 - \rho_{0}}} & {and} & {{\alpha_{1} = \frac{R_{0}\rho_{0}^{2}}{1 - \rho_{1}}},} \end{matrix}$

wherein $\begin{matrix} {R_{0}:=\frac{4\kappa_{i}{{\overset{\sim}{b} - b_{0}}}^{2}}{3d_{i} \times \left( {b_{4} - \overset{\sim}{b}} \right)}} & {and} & {R_{1}:={\frac{4\kappa_{i + 1}{{b_{4} - \overset{\sim}{b}}}^{2}}{3\left( {\overset{\sim}{b} - b_{0}} \right) \times d_{i + 1}}.}} \end{matrix}$


28. A computer program product according to claim 21, wherein the third executable portion determines the interpolant p(t) according to the following: ${{p(t)} = {\sum\limits_{j = 0}^{4}\quad {b_{j}{B_{j}(t)}}}},$

and wherein ${B_{j}(t)}:={\begin{pmatrix} 4 \\ j \end{pmatrix}{{t^{j}\left( {1 - t} \right)}^{({4 - j})}.}}$


29. A computer program product according to claim 21, wherein the second executable portion identifies points b₀, b₁, b₂, b₃ and b₄ such that the quartic interpolant p(t) is sixth-order accurate over the interval (i≦t≦i+1).
 30. A computer program product according to claim 21, wherein the third executable portion determines the quartic interpolant p(t) such that a change in the interpolant p(t) is directly proportional to a change at least one of the provided positions (f_(i), f_(i+1)), directions (d_(i), d_(i+1)) and curvatures (κ_(i), κ_(i+1)). 